#Specify Project Location and File Type

project_location <- '~/Library/CloudStorage/OneDrive-JohnsHopkins/Research Projects/CLIF/CLIF_Projects/ProneProjects/ProningEpi/ProneSevereARF_Output'
file_type <- 'csv'

#Create Sub Folders within Project Folder
# Check if the output directory exists; if not, create it
if (!dir.exists(paste0(project_location, "/summary_analysis"))) {
  dir.create(paste0(project_location, "/summary_analysis"))
}

set.seed(3221984)
df_agg <- read.csv(paste0(project_location, '/summary_output/aggregate_data.csv')) |>
  #Prepare to merge in Hospital Details
  mutate(
    hospital_id=fifelse(
        hospital_id=='HUP', 'penn_hup', hospital_id),
    hospital_id=fifelse(
      hospital_id=='HUPME', 'penn_hupme', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PAH', 'penn_pah', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PPMC', 'penn_ppmc', hospital_id),
    hospital_id=fifelse(
      hospital_id=='UMCOP', 'penn_umcop', hospital_id),
    hospital_id=fifelse(
      hospital_id=='CCH', 'penn_cch', hospital_id)
    )

clif_hospital_description <- read.csv(paste0(project_location, '/clif_hospital_description.csv')) 

#Merge Hospital Data Into DF AGG
df_agg <- df_agg |>
  left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
  #Define Hospital_Size Categorical Variable
  mutate(icu_bed_cat=fcase(
   Number_of_ICU_Beds<20, 0,
   Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
   Number_of_ICU_Beds>=50, 2
  )) |>
 ## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
  mutate(
    icu_bed_cat = fcase(
      Number_of_ICU_Beds < 20,                           0L,      # small
      Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L,      # medium
      Number_of_ICU_Beds >= 50,                           2L,      # large
      default = NA_integer_
    )
  ) %>% 
  ## 2. labelled factor ------------------------------------------------------
  mutate(
    icu_bed_cat = factor(
      icu_bed_cat,
      levels = 0:2,                                    # make sure levels match
      labels = c("small", "medium", "large")
    )
  )
## Joining with `by = join_by(hospital_id)`
#Do this Before Filtering
df_unique_hosp <- df_agg |> distinct(hospital_id, Hospital_Type)
cat('\nHow Many Academic v Community Hospitals are In THese Data')
## 
## How Many Academic v Community Hospitals are In THese Data
table(df_unique_hosp$Hospital_Type)
## 
##  Academic Community 
##        12        25
df_hosptype_summary <- df_agg |>
  group_by(Hospital_Type, study_period) |>
  summarise(
    n_Hospitals=n(),
    n_patients=sum(n_patients, na.rm=T),
    n_proned=sum(observed_prone, na.rm=T),
    percent_proned=round((n_proned/n_patients)*100, 2)
  ) |>
  arrange(study_period, Hospital_Type)
## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.
df_hosptype_summary
## # A tibble: 6 × 6
## # Groups:   Hospital_Type [2]
##   Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
##   <chr>         <chr>              <int>      <int>    <int>          <dbl>
## 1 Academic      COVID                 12       1884      680          36.1 
## 2 Community     COVID                 23       1085      643          59.3 
## 3 Academic      Post-COVID            11       1098      168          15.3 
## 4 Community     Post-COVID            23        498      150          30.1 
## 5 Academic      Pre-COVID             10        982       66           6.72
## 6 Community     Pre-COVID             16        213       30          14.1
#Hospital Size by Period
df_hospsize <- df_agg |>
  group_by(study_period) |>
  summarise(
    small_icus =sum(icu_bed_cat=='small'),
    medium_icus = sum(icu_bed_cat=='medium'),
    large_icus = sum(icu_bed_cat=='large')
  )
cat('\nHospital Size Description')
## 
## Hospital Size Description
df_hospsize
## # A tibble: 3 × 4
##   study_period small_icus medium_icus large_icus
##   <chr>             <int>       <int>      <int>
## 1 COVID                13          10         12
## 2 Post-COVID           12          11         11
## 3 Pre-COVID             8           8         10
#Filter Out Hospitals with Less Than 10 Patients
df_agg <- df_agg |>
  #Filter OUT HOspitals-Periods with < 10 observations
  filter(n_patients>=10)

#Now After Filtering
cat('\nAfer Filtering Out Hospitals with >= 10 Patients \nHow Many Academic v Community Hospitals are In THese Data')
## 
## Afer Filtering Out Hospitals with >= 10 Patients 
## How Many Academic v Community Hospitals are In THese Data
table(df_agg$Hospital_Type)
## 
##  Academic Community 
##        33        50
df_hosptype_summary_post <- df_agg |>
  group_by(Hospital_Type, study_period) |>
  summarise(
    n_Hospitals=n(),
    n_patients=sum(n_patients, na.rm=T),
    n_proned=sum(observed_prone, na.rm=T),
    percent_proned=round((n_proned/n_patients)*100, 2)
  ) |>
  arrange(study_period, Hospital_Type)
## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.
cat('\nProning Counts in Hospitals with >= 10 Patients')
## 
## Proning Counts in Hospitals with >= 10 Patients
df_hosptype_summary
## # A tibble: 6 × 6
## # Groups:   Hospital_Type [2]
##   Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
##   <chr>         <chr>              <int>      <int>    <int>          <dbl>
## 1 Academic      COVID                 12       1884      680          36.1 
## 2 Community     COVID                 23       1085      643          59.3 
## 3 Academic      Post-COVID            11       1098      168          15.3 
## 4 Community     Post-COVID            23        498      150          30.1 
## 5 Academic      Pre-COVID             10        982       66           6.72
## 6 Community     Pre-COVID             16        213       30          14.1
#Hospital Size by Period
df_hospsize <- df_agg |>
  group_by(study_period) |>
  summarise(
    small_icus =sum(icu_bed_cat=='small'),
    medium_icus = sum(icu_bed_cat=='medium'),
    large_icus = sum(icu_bed_cat=='large')
  )
cat('\nHospital Size Description')
## 
## Hospital Size Description
df_hospsize
## # A tibble: 3 × 4
##   study_period small_icus medium_icus large_icus
##   <chr>             <int>       <int>      <int>
## 1 COVID                12          10         12
## 2 Post-COVID            7          11         11
## 3 Pre-COVID             3           7         10
#Academic Size Summary
df_hospsize <- df_agg |>
  group_by(Hospital_Type) |>
  summarise(
    small_icus =sum(icu_bed_cat=='small'),
    medium_icus = sum(icu_bed_cat=='medium'),
    large_icus = sum(icu_bed_cat=='large')
  )
cat('\nHospital Size Description by Hospital Type')
## 
## Hospital Size Description by Hospital Type
df_hospsize
## # A tibble: 2 × 4
##   Hospital_Type small_icus medium_icus large_icus
##   <chr>              <int>       <int>      <int>
## 1 Academic               0           3         30
## 2 Community             22          25          3
#Proning By Hospital Size
df_hospsize <- df_agg |>
  group_by(icu_bed_cat, study_period) |>
  summarise(
    n_Hospitals=n(),
    n_patients=sum(n_patients, na.rm=T),
    n_proned=sum(observed_prone, na.rm=T),
    percent_proned=round((n_proned/n_patients)*100, 2)
  ) |>
  arrange(study_period, icu_bed_cat)
## `summarise()` has grouped output by 'icu_bed_cat'. You can override using the
## `.groups` argument.
cat('\nProning Counts by Hospital Size and Period')
## 
## Proning Counts by Hospital Size and Period
df_hospsize
## # A tibble: 9 × 6
## # Groups:   icu_bed_cat [3]
##   icu_bed_cat study_period n_Hospitals n_patients n_proned percent_proned
##   <fct>       <chr>              <int>      <int>    <int>          <dbl>
## 1 small       COVID                 12        433      300          69.3 
## 2 medium      COVID                 10        623      283          45.4 
## 3 large       COVID                 12       1906      738          38.7 
## 4 small       Post-COVID             7        128       53          41.4 
## 5 medium      Post-COVID            11        335       77          23.0 
## 6 large       Post-COVID            11       1109      179          16.1 
## 7 small       Pre-COVID              3         43        6          14.0 
## 8 medium      Pre-COVID              7        156       23          14.7 
## 9 large       Pre-COVID             10        979       67           6.84
#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned <- df_agg
df_agg_proned <- df_agg_unproned %>%
  mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned <- df_agg_unproned |>
  mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg <- bind_rows(df_agg_proned, df_agg_unproned)

#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
global_agg_byoutcome_period <- read.csv(paste0(project_location, '/summary_output/global_aggregate_byoutcome_period.csv')) |>
  mutate(prone_log_odds = log(prone_rate_adjust / (1 - prone_rate_adjust))) |>
  select(hospital_id, prone_12hour_outcome, prone_log_odds, study_period) |>
  #Prepare to merge in Hospital Details
  mutate(
    hospital_id=fifelse(
        hospital_id=='HUP', 'penn_hup', hospital_id),
    hospital_id=fifelse(
      hospital_id=='HUPME', 'penn_hupme', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PAH', 'penn_pah', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PPMC', 'penn_ppmc', hospital_id),
    hospital_id=fifelse(
      hospital_id=='UMCOP', 'penn_umcop', hospital_id),
    hospital_id=fifelse(
      hospital_id=='CCH', 'penn_cch', hospital_id)
    )

#Merge Back with DF Agg
df_agg <- left_join(df_agg, global_agg_byoutcome_period) |>
  #IF there were no patients proned can use the popensity for unproned patients by filling here
  group_by(hospital_id, study_period) |>
  arrange(hospital_id, study_period, prone_12hour_outcome) |>
  fill(prone_log_odds, .direction = 'down') |>
  ungroup()
## Joining with `by = join_by(hospital_id, study_period, prone_12hour_outcome)`
#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period 
#Also Includes Hospital Type and Hospital Type and Period Interaction
#ICU Capacity is colinear with academic or community so have taken that out of the model

priors <- c(
  #Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in  5-50% range
  prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),  
  #COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
  prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
  #Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
  prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'), 
  #Flat Priors for Hospital Type and Size
  prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
  #Half Normal Weakly Informative but Regularizing Prior for Random Effects
  prior(normal(0,1), class = "sd"))         

#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_full <- brm(
  bf(observed_prone | trials(n_patients) ~ 
       0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type +
     (0 + study_period | hospital_id)), 
  data = df_agg,
  family = binomial,
  prior = priors, 
  cores = 4,
  control = list(adapt_delta = 0.99)
)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
# Examine model fit and summaries
plot(agg_brm_period_full)

summary(agg_brm_period_full)
##  Family: binomial 
##   Links: mu = logit 
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type + (0 + study_period | hospital_id) 
##    Data: df_agg (Number of observations: 166) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36) 
##                                                   Estimate Est.Error l-95% CI
## sd(study_periodCOVID)                                 0.91      0.13     0.70
## sd(study_periodPostMCOVID)                            1.01      0.18     0.72
## sd(study_periodPreMCOVID)                             1.00      0.27     0.58
## cor(study_periodCOVID,study_periodPostMCOVID)         0.77      0.11     0.51
## cor(study_periodCOVID,study_periodPreMCOVID)          0.59      0.21     0.11
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.52      0.23     0.03
##                                                   u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID)                                 1.18 1.00     1198
## sd(study_periodPostMCOVID)                            1.41 1.00     1493
## sd(study_periodPreMCOVID)                             1.61 1.00     1583
## cor(study_periodCOVID,study_periodPostMCOVID)         0.93 1.00     1605
## cor(study_periodCOVID,study_periodPreMCOVID)          0.92 1.00     1599
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.90 1.00     1461
##                                                   Tail_ESS
## sd(study_periodCOVID)                                 2125
## sd(study_periodPostMCOVID)                            2311
## sd(study_periodPreMCOVID)                             2425
## cor(study_periodCOVID,study_periodPostMCOVID)         2070
## cor(study_periodCOVID,study_periodPreMCOVID)          1412
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     1962
## 
## Regression Coefficients:
##                                               Estimate Est.Error l-95% CI
## study_periodCOVID                                -0.11      0.26    -0.61
## study_periodPostMCOVID                           -1.19      0.31    -1.82
## study_periodPreMCOVID                            -2.11      0.34    -2.80
## Hospital_TypeCommunity                            0.85      0.32     0.20
## study_periodPostMCOVID:Hospital_TypeCommunity     0.10      0.31    -0.52
## study_periodPreMCOVID:Hospital_TypeCommunity     -0.16      0.49    -1.11
##                                               u-95% CI Rhat Bulk_ESS Tail_ESS
## study_periodCOVID                                 0.40 1.01      503      873
## study_periodPostMCOVID                           -0.60 1.01      716     1442
## study_periodPreMCOVID                            -1.48 1.00     1242     2148
## Hospital_TypeCommunity                            1.48 1.01      552     1013
## study_periodPostMCOVID:Hospital_TypeCommunity     0.74 1.00     1872     2212
## study_periodPreMCOVID:Hospital_TypeCommunity      0.82 1.00     1797     1795
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
##               3.26 seconds (Total)
## Chain 4: 
## Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 2.033 seconds (Warm-up)
## Chain 3:                2.35 seconds (Sampling)
## Chain 3:                4.383 seconds (Total)
## Chain 3: 
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 2.847 seconds (Warm-up)
## Chain 2:                2.309 seconds (Sampling)
## Chain 2:                5.156 seconds (Total)
## Chain 2:
#Now Generate Odds Ratios for Communtiy vs Academic Proning by Period
full_posterior <- as_draws_df(agg_brm_period_full)  |>
  mutate(pre_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPreMCOVID:Hospital_TypeCommunity`),
         covid = exp(b_Hospital_TypeCommunity),
         post_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPostMCOVID:Hospital_TypeCommunity`))

summary_posterior <- data.frame(
  variable = c(
    'Pre-COVID Community Hospital',
    'COVID Community Hospital', 
    'Post-COVID Community Hospital'
  ),
  median_or = c(median(full_posterior$pre_covid),
              median(full_posterior$covid),
              median(full_posterior$post_covid)),
  conf.low = c(quantile(full_posterior$pre_covid, p=0.025),
               quantile(full_posterior$covid, p=0.025),
               quantile(full_posterior$post_covid, p=0.025)),
  conf.high = c(quantile(full_posterior$pre_covid, p=0.975),
               quantile(full_posterior$covid, p=0.975),
               quantile(full_posterior$post_covid, p=0.975)))
print(summary_posterior)
write.csv(summary_posterior, paste0(project_location, '/summary_output/summary_posterior.csv'))
####Predictions Command from Marginal Effects with 4000 posterior draws
reliability_adjust_full <- predictions(agg_brm_period_full, 
                                  type = 'response', ndraws=4000, re_formula =  ~ (0 + study_period | hospital_id))
reliability_adjust_full <- as_tibble(reliability_adjust_full) |>
  left_join(df_agg %>% 
         select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds', 'prone_log_odds', 'study_period', 'prone_12hour_outcome'), 
         by = c('hospital_id', 'study_period', 'prone_log_odds')) |>
  #Need to Recombine Prone_12hour_outcome Estimates into 1 Per Hospital and Period
  group_by(hospital_id, study_period) |>
  mutate(
    estimate_add=sum(estimate, na.rm=T),
    ci_low_add=sum(conf.low, na.rm=T),
    ci_hi_add=sum(conf.high, na.rm=T),
    n_patients=sum(n_patients, na.rm=T)
  ) |>
  arrange(hospital_id, study_period, -prone_12hour_outcome) |> #THis Arranges so Row Number 1 is the prone_12hour_outcome==1
  filter(row_number()==1) |>
  ungroup() |>
  mutate(adjust_rate=estimate_add/n_patients) |>
  mutate(ci_low=ci_low_add/n_patients) |>
  mutate(ci_hi=ci_hi_add/n_patients) |>
  mutate(period_rank=fcase(
    study_period=='Pre-COVID', 0, 
    study_period=='COVID', 1, 
    study_period=='Post-COVID', 2
  )) |>
  #Orders for Graph
  arrange(period_rank, adjust_rate) |>
  mutate(hospital_rank=row_number()) |>
  mutate(actual_rate=observed_prone/n_patients) |>
  mutate(grouping='Actual Rate') |>
  left_join(df_agg %>% select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds')) |>
  distinct()
## Warning in left_join(as_tibble(reliability_adjust_full), df_agg %>% select("hospital_id", : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 13 of `x` matches multiple rows in `y`.
## ℹ Row 13 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
## Joining with `by = join_by(hospital_id, Number_of_ICU_Beds)`
## Warning in left_join(mutate(mutate(mutate(arrange(mutate(mutate(mutate(mutate(ungroup(filter(arrange(mutate(group_by(left_join(as_tibble(reliability_adjust_full), : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 151 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
##   "many-to-many"` to silence this warning.
break_points <- reliability_adjust_full |>
  mutate(pre_covid_num=fcase(
    row_number()!=1 & study_period=='Pre-COVID' & lead(study_period)=='COVID', hospital_rank),
    post_covid_num=fcase(
      study_period=='COVID' & lead(study_period)=='Post-COVID', hospital_rank
    )) |>
  summarise(
    pre_covid_num=mean(pre_covid_num,na.rm=T)+0.5,
    post_covid_num=mean(post_covid_num, na.rm=T)+0.5)


caterpillar_hospital_type <- ggplot(reliability_adjust_full, 
                              aes(x = hospital_rank,
                                  y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period), color='lightgrey') +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
  geom_point(aes(color=Hospital_Type)) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.1), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals",
    y = "Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals"
  )  +
  geom_vline(xintercept = break_points$pre_covid_num, linetype = 2, linewidth=0.3) +
  annotate("text", x = 7, y = 0.95, label = 'Pre Pandemic', fontface = 'bold') +
  geom_vline(xintercept = break_points$post_covid_num, linetype = 2, linewidth=0.3) +
  annotate("text", x = 35.5, y = 0.95, label = 'Pandemic', fontface = 'bold') +
  annotate("text", x = 68, y = 0.95, label = 'Post Pandemic', fontface = 'bold') +
  theme_classic() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank(),
        axis.line.x = element_line(color = "black", linewidth = 0.5)) 
caterpillar_hospital_type

ggsave('caterpillar_hospital_byperiod_fully_adjusted.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
  height = 6,       # adjust height for PowerPoint-friendly proportions
  units = "in")    # units set to inches for precise control)
  
#With Actual Rate
#Add the Actual Value
with_rate <- caterpillar_hospital_type +
  geom_point(aes(y=observed_prone/n_patients), color = 'maroon', alpha = 0.25)
with_rate

ggsave('caterpillar_fully_adjusted.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
  height = 6,       # adjust height for PowerPoint-friendly proportions
  units = "in")    # units set to inches for precise control)
#Stack Caterpillar Plots Vertically
vertical_data <- reliability_adjust_full |>
  group_by(study_period) |>
  arrange(study_period, adjust_rate) |>
  mutate(hospital_rank=row_number()) |>
  ungroup()

vertical <-ggplot(vertical_data, 
                              aes(x = hospital_rank,
                                  y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period), color='lightgrey') +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
  geom_point(aes(color=Hospital_Type)) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.2), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals\nRanked by Proning Use",
    y = "% Eligible Patients Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals"
  )  +
  theme_minimal() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank(),
        axis.line.x = element_line(color = "black", linewidth = 0.25),
        strip.background = element_rect(fill = "grey90", color = NA),
    strip.text       = element_text(color = "black"),
    legend.position = c(0.1, 0.90),
    legend.background = element_rect(fill = "white", color = "grey"),
    legend.title = element_text(size = 9)) +
    facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')), 
               ncol=1,
               scales='free_x')
## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ggsave('vertical_caterpillar.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
  height = 6,       # adjust height for PowerPoint-friendly proportions
  units = "in")    # units set to inches for precise control)

verticald <-ggplot(vertical_data, 
                  aes(x = hospital_rank,
                      y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period), color='lightgrey') +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
  geom_point(aes(color=Hospital_Type)) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.2), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals\nRanked by Proning Use",
    y = "% Eligible Patients Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals"
  )  +
  theme_minimal() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank(),
        axis.line.x = element_line(color = "black", linewidth = 0.25),
        strip.background = element_rect(fill = "grey90", color = NA),
        strip.text       = element_text(color = "black"),
        legend.position = c(0.75, 0.85),
        legend.background = element_rect(fill = "white", color = "grey"),
        legend.title = element_text(size = 12),
        legend.key.height   = unit(1.5, "lines"),
        legend.key.width = unit(2.2, 'lines')) +
  facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')), 
             ncol=1)

ggsave('vertical_caterpillar_option_D.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
       height = 6,       # adjust height for PowerPoint-friendly proportions
       units = "in")    # units set to inches for precise control)

stacked_data <- vertical_data |>
  group_by(study_period) |>
  arrange(study_period, adjust_rate) |>
  mutate(hospital_rank=if_else(
    study_period=='COVID', row_number(), NA_real_
  )) |>
  group_by(hospital_id) |>
  fill(hospital_rank, .direction = 'updown')

stacked <- ggplot(stacked_data, 
                  aes(x = hospital_rank,
                      y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period, color = study_period)) +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= study_period), alpha = 0.4) +
  geom_point(aes(color=study_period, shape=Hospital_Type), size = 2.5) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.2), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals\nRanked by COVID Period Proning Use",
    y = "% Eligible Patients Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals"
  )  +
  theme_minimal() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank())

ggsave('stacked_caterpillar.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
       height = 6,       # adjust height for PowerPoint-friendly proportions
       units = "in")    # units set to inches for precise control)
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
stacked2 <- ggplot(stacked_data, 
                              aes(x = hospital_rank,
                                  y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period, linetype=study_period),  color = 'black') +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) +
  geom_point(aes(color=Hospital_Type, shape=study_period), size = 2) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.2), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Demuth") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals\nRanked by COVID Period Proning Use",
    y = "% Eligible Patients Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals"
  )  +
  theme_minimal() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank(),
        legend.position = c(0.15, 0.75),
        legend.background = element_rect(fill = "white", color = "grey"))

ggsave('stacked_caterpillar2.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
       height = 6,       # adjust height for PowerPoint-friendly proportions
       units = "in")    # units set to inches for precise control)
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
stacked_data <- stacked_data |>
  group_by(study_period) |>
  mutate(period_n=sum(n_patients)) |>
  ungroup()

stacked3 <- ggplot(stacked_data, 
                   aes(x = hospital_rank,
                       y = adjust_rate)) +
  geom_line(aes(x = hospital_rank,
                y = adjust_rate, group=study_period, alpha=period_n),  color = 'black') +
  geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) + scale_alpha(guide="none")+
  geom_point(aes(color=Hospital_Type, shape=study_period), size = 2.25) +  # Add points for center effects
  scale_x_continuous(
    breaks = seq(1, nrow(reliability_adjust_full), by=1), 
    labels = NULL  # Optionally use hospital IDs as labels
  ) +
  scale_y_continuous(breaks=seq(0,1.0, by=0.2), 
                     labels = scales::percent, 
                     limits = c(0,1.0)) +
  MetBrewer::scale_color_met_d("Egypt") + #MetBrewer has many pallettes to choose from
  labs(
    x = "Hospitals\nRanked by Pandemic Period Proning Use",
    y = "% Eligible Patients Proned Within 12 Hours",
    caption = "Error bars represent 95% credible intervals", 
    shape="Study Period", 
    color= "Hospital Type"
  )  +
  scale_shape_manual(values= c("circle", "triangle", "square"), labels=c("Pandemic", "Post-pandemic", "Pre-pandemic"))+
  theme_classic() +
  theme(panel.grid.major.x=element_blank(),
        panel.grid.minor.x=element_blank(),
        panel.grid.minor.y=element_blank(),
        legend.position = c(0.15, 0.75),
        legend.background = element_rect(fill = "white", color = "grey"))

ggsave('stacked_caterpillar3.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'),
       width = 12,       # adjust width to make it wide
       height = 6,       # adjust height for PowerPoint-friendly proportions
       units = "in")    # units set to inches for precise control)
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
posterior <- as_draws_df(agg_brm_period_full)

# Extract the standard deviations (SDs) for each study_period random effect
sd_pre <- posterior$`sd_hospital_id__study_periodPreMCOVID`
sd_covid <- posterior$`sd_hospital_id__study_periodCOVID`
sd_post <- posterior$`sd_hospital_id__study_periodPostMCOVID`

#3 Calculate MOR for Each Period
# MOR calculation function
calc_mor <- function(sd_samples) {
  exp(sqrt(2) * sd_samples * qnorm(0.75))  # qnorm(0.75) ~ 0.674
}

# Compute MORs
mor_pre <- calc_mor(sd_pre)
mor_covid <- calc_mor(sd_covid)
mor_post <- calc_mor(sd_post)

# Differences Using The Latter Periods as Reference
diff_covid_pre <-  mor_covid - mor_pre
diff_post_pre <- mor_post - mor_pre
diff_post_covid <- mor_post - mor_covid

# Posterior probabilities that the Difference in mOR is > 0 (Bayesian 'p-values')
prop_covid_gt_pre <- mean(diff_covid_pre  > 0)
prop_post_gt_pre  <- mean(diff_post_pre > 0)
prop_post_gt_covid <- mean(diff_post_covid > 0)

# Credible intervals
cr_covid_pre <- quantile(diff_covid_pre, c(0.025, 0.5, 0.975))
cr_post_pre <- quantile(diff_post_pre, c(0.025, 0.5, 0.975))
cr_post_covid <- quantile(diff_post_covid, c(0.025, 0.5, 0.975))

# Summarize
mor_summary <- data.frame(
  Period = c("Pre-MCOVID", "COVID", "Post-MCOVID"),
  MOR_Mean = c(mean(mor_pre), mean(mor_covid), mean(mor_post)),
  MOR_Median = c(median(mor_pre), median(mor_covid), median(mor_post)),
  MOR_2.5 = c(quantile(mor_pre, 0.025), quantile(mor_covid, 0.025), quantile(mor_post, 0.025)),
  MOR_97.5 = c(quantile(mor_pre, 0.975), quantile(mor_covid, 0.975), quantile(mor_post, 0.975))
)

print(mor_summary)
##        Period MOR_Mean MOR_Median  MOR_2.5 MOR_97.5
## 1  Pre-MCOVID 2.690205   2.517063 1.745838 4.657248
## 2       COVID 2.396583   2.354657 1.944955 3.090758
## 3 Post-MCOVID 2.670977   2.587785 1.990022 3.836999
#Summarize the Differences in Median Odds Ratios by Period
diff_mor <- data.frame(
  contrast = c('COVID-Pre', 'Post-COVID', 'Post-Pre'),
  MOR_diff_mean = c(mean(diff_covid_pre), mean(diff_post_covid), mean(diff_post_pre)),
  MOR_diff_median = c(median(diff_covid_pre), median(diff_post_covid), median(diff_post_pre)),
  MOR_diff_2.5 = c(quantile(diff_covid_pre, 0.025), quantile(diff_post_covid, 0.025), quantile(diff_post_pre, 0.025)),
  MOR_diff_97.5 = c(quantile(diff_covid_pre, 0.975), quantile(diff_post_covid, 0.975), quantile(diff_post_pre, 0.975)),
  prob_diff_gt_0 = c(prop_covid_gt_pre, prop_post_gt_covid, prop_post_gt_pre)
)
print(diff_mor)
##     contrast MOR_diff_mean MOR_diff_median MOR_diff_2.5 MOR_diff_97.5
## 1  COVID-Pre   -0.29362135     -0.15911771   -2.2612521     0.8422212
## 2 Post-COVID    0.27439383      0.22563714   -0.5031843     1.3230072
## 3   Post-Pre   -0.01922753      0.07103335   -2.0632630     1.4426286
##   prob_diff_gt_0
## 1        0.39675
## 2        0.73075
## 3        0.55000
#Visualize
mor_df <- data.frame(
  MOR = c(mor_covid, mor_pre, mor_post),
  Period = factor(rep(c("COVID", "Pre-MCOVID", "Post-MCOVID"), each = length(mor_covid)))
)

ggplot(mor_df, aes(x = MOR, fill = Period)) +
  geom_density(alpha = 0.5) +
  labs(
    title = "Period-Specific Median Odds Ratios (MOR  \n (Fully Adjusted))",
    x = "Median Odds Ratio (MOR)",
    y = "Density"
  ) +
  scale_x_continuous(seq(1,20, by=1), limits =c(1,20),
                     name='Distribution of Median Odds Ratos by Study Period') +
  theme_minimal()

ggsave('mor_distribution_fully.pdf', 
       path = paste0(project_location, '/summary_output/graphs/'))
## Saving 7 x 5 in image
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period 
#Adjusted for COVID Status

#First Prepare Dataset
#Open Aggregate Data with COVID
df_agg_covid <- read.csv(paste0(project_location, '/summary_output/aggregate_data_wcovid.csv')) |>
  #Prepare to merge in Hospital Details
  mutate(
    hospital_id=fifelse(
        hospital_id=='HUP', 'penn_hup', hospital_id),
    hospital_id=fifelse(
      hospital_id=='HUPME', 'penn_hupme', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PAH', 'penn_pah', hospital_id),
    hospital_id=fifelse(
      hospital_id=='PPMC', 'penn_ppmc', hospital_id),
    hospital_id=fifelse(
      hospital_id=='UMCOP', 'penn_umcop', hospital_id),
    hospital_id=fifelse(
      hospital_id=='CCH', 'penn_cch', hospital_id)
    )

#Merge Hospital Data Into DF AGG
df_agg_covid <- df_agg_covid |>
  left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
  #Define Hospital_Size Categorical Variable
  mutate(icu_bed_cat=fcase(
   Number_of_ICU_Beds<20, 0,
   Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
   Number_of_ICU_Beds>=50, 2
  )) |>
 ## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
  mutate(
    icu_bed_cat = fcase(
      Number_of_ICU_Beds < 20,                           0L,      # small
      Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L,      # medium
      Number_of_ICU_Beds >= 50,                           2L,      # large
      default = NA_integer_
    )
  ) %>% 
  ## 2. labelled factor ------------------------------------------------------
  mutate(
    icu_bed_cat = factor(
      icu_bed_cat,
      levels = 0:2,                                    # make sure levels match
      labels = c("small", "medium", "large")
    )
  )
## Joining with `by = join_by(hospital_id)`
#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned_covid <- df_agg_covid
df_agg_proned_covid <- df_agg_unproned_covid %>%
  mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned_covid <- df_agg_unproned_covid |>
  mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg_covid <- bind_rows(df_agg_proned_covid, df_agg_unproned_covid) |>
  #Only have the propensity at the period and outcome level so need to add period back here
  mutate(study_period=fcase(
    period_sarscov2=='COVID_SarsCov2Neg-Unk', 'COVID',
    period_sarscov2=='COVID_SarsCov2Pos', 'COVID',
    period_sarscov2=='Post-COVID_SarsCov2Neg-Unk', 'Post-COVID',
    period_sarscov2=='Post-COVID_SarsCov2Pos', 'Post-COVID',
    period_sarscov2=='Pre-COVID', 'Pre-COVID'
  ))


#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
#Merge Back with DF Agg COVID
df_agg_covid <- left_join(df_agg_covid, global_agg_byoutcome_period) |>
  #IF there were no patients proned can use the popensity for unproned patients by filling here
  group_by(hospital_id, period_sarscov2) |>
  arrange(hospital_id, period_sarscov2, prone_12hour_outcome) |>
  fill(prone_log_odds, .direction = 'down') |>
  ungroup()
## Joining with `by = join_by(hospital_id, prone_12hour_outcome, study_period)`
set.seed(32284)

priors <- c(
  #Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in  5-50% rnage
  prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PreMCOVID'),  
  #COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
  prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2Pos'),
  prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2NegMUnk'),
  #Post-COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
  prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2Pos'), 
  prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2NegMUnk'), 
  #Flat Priors for Hospital Type and Size
  prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
  #Half Normal Weakly Informative but Regularizing Prior for Random Effects
  prior(normal(0,1), class = "sd"))         

#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_covid <- brm(
  bf(observed_prone | trials(n_patients) ~ 
       0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type +
       (0 + period_sarscov2 | hospital_id)), 
  data = df_agg_covid,
  family = binomial,
  prior = priors, 
  cores = 4,
  control = list(adapt_delta = 0.99)
)
## Warning: Rows containing NAs were excluded from the model.
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
plot(agg_brm_period_covid)

summary(agg_brm_period_covid)
##  Family: binomial 
##   Links: mu = logit 
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type + (0 + period_sarscov2 | hospital_id) 
##    Data: df_agg_covid (Number of observations: 318) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 37) 
##                                                                                      Estimate
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                 0.94
## sd(period_sarscov2COVID_SarsCov2Pos)                                                     0.91
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                            0.97
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                0.61
## sd(period_sarscov2PreMCOVID)                                                             0.96
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)               0.86
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      0.79
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)          0.68
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)          0.56
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)              0.52
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)     0.62
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                       0.48
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                           0.57
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                  0.42
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                      0.39
##                                                                                      Est.Error
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                  0.15
## sd(period_sarscov2COVID_SarsCov2Pos)                                                      0.12
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                             0.16
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                 0.28
## sd(period_sarscov2PreMCOVID)                                                              0.25
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)                0.09
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)       0.12
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)           0.13
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)           0.27
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)               0.26
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)      0.26
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                        0.21
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                            0.19
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                   0.22
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                       0.30
##                                                                                      l-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                 0.67
## sd(period_sarscov2COVID_SarsCov2Pos)                                                     0.69
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                            0.70
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                0.10
## sd(period_sarscov2PreMCOVID)                                                             0.57
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)               0.64
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      0.51
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)          0.37
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)         -0.09
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)             -0.10
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)    -0.05
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                       0.01
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                           0.14
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                 -0.05
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                     -0.26
##                                                                                      u-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                 1.28
## sd(period_sarscov2COVID_SarsCov2Pos)                                                     1.18
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                            1.31
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                1.23
## sd(period_sarscov2PreMCOVID)                                                             1.54
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)               0.97
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      0.95
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)          0.88
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)          0.92
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)              0.90
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)     0.94
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                       0.84
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                           0.88
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                  0.80
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                      0.87
##                                                                                      Rhat
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                             1.00
## sd(period_sarscov2COVID_SarsCov2Pos)                                                 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                        1.00
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                            1.00
## sd(period_sarscov2PreMCOVID)                                                         1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)           1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)  1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)      1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)          1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                   1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                       1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)              1.00
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                  1.00
##                                                                                      Bulk_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                 1395
## sd(period_sarscov2COVID_SarsCov2Pos)                                                     1947
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                            2339
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                2033
## sd(period_sarscov2PreMCOVID)                                                             1923
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)               1359
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      1448
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)          2034
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)          3045
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)              3958
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)     3342
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                       1887
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                           2653
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                  2735
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                      1198
##                                                                                      Tail_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk)                                                 2121
## sd(period_sarscov2COVID_SarsCov2Pos)                                                     2828
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk)                                            2864
## sd(period_sarscov2PostMCOVID_SarsCov2Pos)                                                1504
## sd(period_sarscov2PreMCOVID)                                                             2689
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos)               1861
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk)      2282
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk)          3171
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)          2448
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos)              2573
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos)     2498
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                       2228
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID)                           2631
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID)                  3061
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID)                      2439
## 
## Regression Coefficients:
##                                           Estimate Est.Error l-95% CI u-95% CI
## period_sarscov2COVID_SarsCov2NegMUnk         -0.99      0.22    -1.43    -0.56
## period_sarscov2COVID_SarsCov2Pos              0.46      0.22     0.01     0.88
## period_sarscov2PostMCOVID_SarsCov2NegMUnk    -1.25      0.24    -1.71    -0.78
## period_sarscov2PostMCOVID_SarsCov2Pos        -0.27      0.27    -0.83     0.25
## period_sarscov2PreMCOVID                     -2.29      0.29    -2.89    -1.78
## Hospital_TypeCommunity                        0.84      0.25     0.34     1.33
##                                           Rhat Bulk_ESS Tail_ESS
## period_sarscov2COVID_SarsCov2NegMUnk      1.01      753     1254
## period_sarscov2COVID_SarsCov2Pos          1.01      846     1578
## period_sarscov2PostMCOVID_SarsCov2NegMUnk 1.00      975     1800
## period_sarscov2PostMCOVID_SarsCov2Pos     1.00     1768     2496
## period_sarscov2PreMCOVID                  1.00     1617     2477
## Hospital_TypeCommunity                    1.01      973     1833
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
#Compare Differences in COVID vs Not COVID in COVID and POST COVID Period
post_covid <- as_draws_df(agg_brm_period_covid)

# Compute OR for the contrast
delta_covid_period <- post_covid$b_period_sarscov2COVID_SarsCov2Pos - post_covid$b_period_sarscov2COVID_SarsCov2NegMUnk
or <- exp(delta_covid_period)

# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period \n', posterior_summary)
## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period 
##  4.2533 3.20544 5.561699
# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2Pos - post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk
or <- exp(delta_post_covid_period)

# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period \n', posterior_summary)
## 
## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period 
##  2.666524 1.570643 4.490486
#COVID Negative in COVID Period vs Pre-COVID
# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk - post_covid$b_period_sarscov2PreMCOVID
or <- exp(delta_post_covid_period)

# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID \n', posterior_summary)
## 
## Odds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID 
##  2.812252 1.649069 5.135836
#ICU bed capacity Not Included as an additional covariable (in addition to hospital type becaues of co-linearity)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period 
#ICU Capacity is colinear with academic or community so have taken that out of the model

priors <- c(
  #Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in  5-50% range
  prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),  
  #COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
  prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
  #Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
  prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'), 
  #Flat Priors for Hospital Size
  prior(normal(0,1), class = "b", coef = 'icu_bed_catlarge'),
  prior(normal(0,1), class = "b", coef = 'icu_bed_catmedium'),
  #Half Normal Weakly Informative but Regularizing Prior for Random Effects
  prior(normal(0,1), class = "sd"))         

#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_size <- brm(
  bf(observed_prone | trials(n_patients) ~ 
       0 + offset(prone_log_odds) + study_period + icu_bed_cat +
     (0 + study_period | hospital_id)), 
  data = df_agg,
  family = binomial,
  prior = priors, 
  cores = 4,
  control = list(adapt_delta = 0.99)
)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
# Examine model fit and summaries
plot(agg_brm_period_size)

summary(agg_brm_period_size)
##  Family: binomial 
##   Links: mu = logit 
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + icu_bed_cat + (0 + study_period | hospital_id) 
##    Data: df_agg (Number of observations: 166) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36) 
##                                                   Estimate Est.Error l-95% CI
## sd(study_periodCOVID)                                 0.93      0.13     0.71
## sd(study_periodPostMCOVID)                            1.03      0.18     0.72
## sd(study_periodPreMCOVID)                             1.02      0.26     0.60
## cor(study_periodCOVID,study_periodPostMCOVID)         0.80      0.10     0.54
## cor(study_periodCOVID,study_periodPreMCOVID)          0.63      0.19     0.18
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.54      0.22     0.06
##                                                   u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID)                                 1.24 1.00     1224
## sd(study_periodPostMCOVID)                            1.42 1.00     1581
## sd(study_periodPreMCOVID)                             1.61 1.00     1858
## cor(study_periodCOVID,study_periodPostMCOVID)         0.94 1.00     1732
## cor(study_periodCOVID,study_periodPreMCOVID)          0.93 1.00     1928
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.91 1.00     1403
##                                                   Tail_ESS
## sd(study_periodCOVID)                                 1816
## sd(study_periodPostMCOVID)                            2556
## sd(study_periodPreMCOVID)                             2671
## cor(study_periodCOVID,study_periodPostMCOVID)         2488
## cor(study_periodCOVID,study_periodPreMCOVID)          2396
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     2180
## 
## Regression Coefficients:
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID          0.67      0.25     0.16     1.16 1.00     1174
## study_periodPostMCOVID    -0.35      0.28    -0.90     0.19 1.00     1309
## study_periodPreMCOVID     -1.42      0.33    -2.11    -0.78 1.00     1223
## icu_bed_catmedium         -0.32      0.36    -1.00     0.39 1.00     1232
## icu_bed_catlarge          -0.55      0.33    -1.19     0.08 1.00      977
##                        Tail_ESS
## study_periodCOVID          1917
## study_periodPostMCOVID     2112
## study_periodPreMCOVID      2316
## icu_bed_catmedium          2041
## icu_bed_catlarge           1667
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## ing)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.302 seconds (Warm-up)
## Chain 1:                2.729 seconds (Sampling)
## Chain 1:                5.031 seconds (Total)
## Chain 1: 
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.822 seconds (Warm-up)
## Chain 4:                2.316 seconds (Sampling)
## Chain 4:                5.138 seconds (Total)
## Chain 4: 
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 3.078 seconds (Warm-up)
## Chain 3:                3.011 seconds (Sampling)
## Chain 3:                6.089 seconds (Total)
## Chain 3: 
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 3.266 seconds (Warm-up)
## Chain 2:                3.066 seconds (Sampling)
## Chain 2:                6.332 seconds (Total)
## Chain 2:
#Now Generate Odds Ratios for Hospital Size Contrasts
full_posterior <- as_draws_df(agg_brm_period_size)  |>
  mutate(medium_v_small = exp(b_icu_bed_catmedium ),
         large_v_small = exp(b_icu_bed_catlarge))

summary_posterior_size <- data.frame(
  variable = c(
    'Medium vs Small ICU Capacity',
    'Large vs Small ICU Capacity'
  ),
  median_or = c(median(full_posterior$medium_v_small),
              median(full_posterior$large_v_small)),
  conf.low = c(quantile(full_posterior$medium_v_small, p=0.025),
               quantile(full_posterior$large_v_small, p=0.025)),
  conf.high = c(quantile(full_posterior$medium_v_small, p=0.975),
               quantile(full_posterior$large_v_small, p=0.975)))
print(summary_posterior_size)
write.csv(summary_posterior_size, paste0(project_location, '/summary_output/summary_posterior_size_analysis.csv'))
#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period 

set.seed(32284)

priors <- c(
  #Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in  5-50% rnage
  prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),  
  #COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
  prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
  #Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
  prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'), 
  #Half Normal Weakly Informative but Regularizing Prior for Random Effects
  prior(normal(0,1), class = "sd"))         

#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_strat <- brm(
  bf(observed_prone | trials(n_patients) ~ 
       0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)), 
  data = df_agg,
  family = binomial,
  prior = priors, 
  cores = 4,
  control = list(adapt_delta = 0.99)
)
## Compiling Stan program...
## Trying to compile a simple C file
## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
##   679 | #include <cmath>
##       |          ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1
## Start sampling
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
# Examine model fit and summaries
plot(agg_brm_period_strat)

summary(agg_brm_period_strat)
##  Family: binomial 
##   Links: mu = logit 
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id) 
##    Data: df_agg (Number of observations: 166) 
##   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup draws = 4000
## 
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36) 
##                                                   Estimate Est.Error l-95% CI
## sd(study_periodCOVID)                                 0.98      0.13     0.77
## sd(study_periodPostMCOVID)                            1.09      0.18     0.79
## sd(study_periodPreMCOVID)                             1.03      0.25     0.63
## cor(study_periodCOVID,study_periodPostMCOVID)         0.82      0.09     0.61
## cor(study_periodCOVID,study_periodPreMCOVID)          0.65      0.19     0.22
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.57      0.22     0.11
##                                                   u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID)                                 1.27 1.00      790
## sd(study_periodPostMCOVID)                            1.50 1.00     1292
## sd(study_periodPreMCOVID)                             1.59 1.00     1646
## cor(study_periodCOVID,study_periodPostMCOVID)         0.95 1.00     1480
## cor(study_periodCOVID,study_periodPreMCOVID)          0.94 1.00     1424
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     0.92 1.00     1306
##                                                   Tail_ESS
## sd(study_periodCOVID)                                 1734
## sd(study_periodPostMCOVID)                            2526
## sd(study_periodPreMCOVID)                             2829
## cor(study_periodCOVID,study_periodPostMCOVID)         1909
## cor(study_periodCOVID,study_periodPreMCOVID)          1739
## cor(study_periodPostMCOVID,study_periodPreMCOVID)     1794
## 
## Regression Coefficients:
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID          0.39      0.17     0.05     0.72 1.02      362
## study_periodPostMCOVID    -0.64      0.20    -1.06    -0.26 1.01      688
## study_periodPreMCOVID     -1.72      0.26    -2.24    -1.23 1.01     1103
##                        Tail_ESS
## study_periodCOVID          1009
## study_periodPostMCOVID     1469
## study_periodPreMCOVID      2331
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## ration: 1000 / 2000 [ 50%]  (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 1.733 seconds (Warm-up)
## Chain 1:                1.942 seconds (Sampling)
## Chain 1:                3.675 seconds (Total)
## Chain 1: 
## Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.826 seconds (Warm-up)
## Chain 2:                1.936 seconds (Sampling)
## Chain 2:                3.762 seconds (Total)
## Chain 2: 
## Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 1.791 seconds (Warm-up)
## Chain 3:                2.196 seconds (Sampling)
## Chain 3:                3.987 seconds (Total)
## Chain 3: 
## Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.618 seconds (Warm-up)
## Chain 4:                1.928 seconds (Sampling)
## Chain 4:                4.546 seconds (Total)
## Chain 4:
#Posterior Draws
sensitivity_posterior <- as_draws_df(agg_brm_period_strat) |>
  mutate(covid_v_pre_covid=exp(b_study_periodCOVID-b_study_periodPreMCOVID),
         covid_v_post_covid=exp(b_study_periodCOVID-b_study_periodPostMCOVID),
         pre_covid_v_covid=exp(b_study_periodPreMCOVID-b_study_periodCOVID),
         post_covid_v_covid=exp(b_study_periodPostMCOVID-b_study_periodCOVID),
         post_covid_v_precovid=exp(b_study_periodPostMCOVID-b_study_periodPreMCOVID))

sensitivity_table <- data.frame(
  variable = c('COVID_v_Pre-COVID',
               'COVID_v_Post-COVID',
               'Pre-COVID_v_COVID',
               'Post-COVID_v_COVID',
               'Post-COVID_v_Pre-COVID'),
  or = c(median(sensitivity_posterior$covid_v_pre_covid),
         median(sensitivity_posterior$covid_v_post_covid),
         median(sensitivity_posterior$pre_covid_v_covid),
         median(sensitivity_posterior$post_covid_v_covid),
         median(sensitivity_posterior$post_covid_v_precovid)),
  conf.low=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.025),
         quantile(sensitivity_posterior$covid_v_post_covid, p=0.025),
         quantile(sensitivity_posterior$pre_covid_v_covid, p=0.025),
         quantile(sensitivity_posterior$post_covid_v_covid, p=0.025),
         quantile(sensitivity_posterior$post_covid_v_precovid, p=0.025)),
  conf.high=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.975),
         quantile(sensitivity_posterior$covid_v_post_covid, p=0.975),
         quantile(sensitivity_posterior$pre_covid_v_covid, p=0.975),
         quantile(sensitivity_posterior$post_covid_v_covid, p=0.975),
         quantile(sensitivity_posterior$post_covid_v_precovid, p=0.975)))
sensitivity_table
write.csv(sensitivity_table, paste0(project_location, '/summary_output/sensitivity_table.csv'))